/**
 Configuration.cxx
 \file Implementation of class Configuration.
 
 \author Thomas Burton
 \date 11/17/11
 */

#include "tppmc/Configuration.h"

#include <cmath>
#include <cstdlib>
#include <sstream>
#include <stdexcept>

#include <getopt.h>

#include <boost/foreach.hpp>

#include "tppmc/Enumerations.h"

namespace tppmc {

const double PYTHIA_MIN_PARTON_PT(2.5);

enum MachineType{RHIC_200,RHIC_500,FNAL_E704_PROTON,FNAL_E704_PBAR};
enum GluonSiversType{NO_SIVERS,AS_UP_QUARK,POS_LIMIT,AS_UP_QUARK_POS_BND};
enum PythiaEventType{QCD2x2,GENJPSI};

Configuration::Configuration(int argc, char* argv[])
: mPytune(0)
, mLhapdf(-1, -1) {
   fftype  = 4;  
   order   = LO; 
   numEvts = 0; 
   machine_type = RHIC_200; 
   pyth_evt_type = QCD2x2;
   quark_transv_soffer = 0; 
   quark_siv_posbnd = 0; 
   etamin = -999.;
   etamax = 999.;
   mean_jt2 = 0.25;
   min_hadron_pt = 0.;
   fav_pos_bnd = 0; 
   unfav_pos_bnd = 0; 
   helpflag = 0; 
   sivers_type = 0;
   pyth_parton_pt = PYTHIA_MIN_PARTON_PT;
   primordial_kt = std::sqrt(0.25);
   // Default process selection for QCD 2 --> 2 scattering
   int processes[7] = {11, 12, 13, 28, 53, 68, 96};
   mProcesses = std::set<int>(processes, processes + 7);
   tparticleout = 0;
   if(argv) {
      Set(argc, argv);
   } // if
}

Configuration::~Configuration() {
}

int Configuration::Set(int argc, char** argv) {
   /* no_argument, required_argument, optind, option are all defined
    in getopt.h.
    struct option {
    char* name;
    int has_arg;
    int* flag;
    int val;
    };
    */
   option long_options[] = {
      {"help",                   no_argument, &helpflag, 1 },
      {"quark_transv_soffer",    no_argument, &quark_transv_soffer, 1 },
      {"tppmc_frag_function",    required_argument, NULL, 'g'},
      {"ffunc",                  required_argument, NULL, 'f'},
      {"output",                 required_argument, NULL, 'o'},
      {"order",                  required_argument, NULL, 'r'},
      {"events",                 required_argument, NULL, 'e'},
      {"rhic200",                no_argument, &machine_type, RHIC_200 },
      {"rhic500",                no_argument, &machine_type, RHIC_500 },
      {"fnalftp",                no_argument, &machine_type, FNAL_E704_PROTON },
      {"fnalftpb",               no_argument, &machine_type, FNAL_E704_PBAR },
      {"qcd2x2",                 no_argument, &pyth_evt_type, QCD2x2},
      {"processes",              required_argument, NULL, 'P'},
      {"min_hard_pt",            required_argument, NULL, 'm'},
      {"jpsi",                   no_argument, &pyth_evt_type, GENJPSI},
      {"pytune",                 required_argument, NULL, 'p'},
      {"qa",                     required_argument, NULL, 'q'},
      {"lhapdf",                 required_argument, NULL, 'l'},
      {"sivers",                 required_argument, NULL, 's'},
      {"eta",                    required_argument, NULL, 'E'},
      {"meanjt2",                required_argument, NULL, 'j'},
      {"min_hadron_pt",          required_argument, NULL, 't'},
      {"primordial_kt",          required_argument, NULL, 'k'},
      {"hadron",                 required_argument, NULL, 'h'},
      {"tparticleout",           required_argument, NULL, 'n'},
      {"transv",                 required_argument, NULL, 'T'},
      {NULL, 0, NULL, 0}
   };
   while(true){
      int option_index(0);
      int c = getopt_long(argc, argv,"g:f:o:r:e:p:P:l:q:s:E:j:t:m:k:h:n:T:",
                          long_options, &option_index);
      if(c == -1) {
         break;
      } // if
      std::stringstream ss; // For string --> number conversions
      int integer(0);         // Integer temporary for use in switch statement
      // Process short options
      switch(c) {
         case 0:
            if(long_options[option_index].flag != 0) {
               break;
            } // if
            printf("option %s", long_options[option_index].name);
            if(optarg) {
               printf(" with arg %s", optarg);
            } // if
            printf("\n");
            break;
         case 'f':
            fftype = strtol(optarg, NULL, 10);
            break;
         case 'o':
            outfile = optarg; 
            break;
         case 'r':
            order = strtol(optarg, NULL, 10);
            break;
         case 'e':
            numEvts = strtol(optarg, NULL, 10);
            break;
         case 'n':
            tparticleout = strtol(optarg, NULL, 10);
            break;
         case 'k':
            primordial_kt = strtof(optarg, NULL);
            break;
         case 'p':
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> mPytune;
            break;
         case 'P':
            ss.str("");
            ss.clear();
            ss << optarg;
            mProcesses.clear();
            while(ss.good()) {
               ss >> integer;
               mProcesses.insert(integer);
            } // while
            break;
         case 'l':
#ifdef USING_LHAPDF
            // Input format is "set,member"
            std::cout << optarg << std::endl;
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> mLhapdf.first >> mLhapdf.second;
#else
            std::cerr <<
            "LHAPDF functionality is not supported." <<
            "\nRun configure --help for how to compile with LHAPDF support."
            << std::endl;
#endif
            break;
         case 'q':
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> mQaFiguresFileName;
            break;
         case 's':
            sivers_type = strtol(optarg, NULL, 10);
            break;
         case 'E':
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> etamin >> etamax;
            break;
         case 'j':
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> mean_jt2;
            break;
         case 't':
            std::cout << optarg << std::endl;
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> min_hadron_pt;
            break;
         case 'm':
            pyth_parton_pt = strtof(optarg, NULL);
            break; 
         case 'h':
            // Push the next hadron type to the set.
            ss.str("");
            ss.clear();
            ss << optarg;
            ss >> integer;
            if(not ss.fail()) {
               mHadrons.insert(integer);
            } // if
            break;
         case 'T':
            transversity_type = strtol(optarg, NULL, 10);
            std::cout << "Selected transversity function " << transversity_type << std::endl;
            break;
         default:
            abort();
      } // switch
   } // while
   if(optind < argc or 1 == helpflag) {
      PrintUsage(std::cout);
   } // if
   if((0 == numEvts or outfile.empty()) and mQaFiguresFileName.empty()) {
      std::cout <<
      "Must specify at least an output file (-o) and number of events (-e) to run." << std::endl;
      exit(-1); 
   } // if
   if(not Verify()) {
      throw std::runtime_error("Invalid Configuration");
   } // if
   return optind;
}

void Configuration::PrintUsage(std::ostream& o) {
   o << std::endl; 
   o << "usage: tppmc -e [numevents] -o [output file]\n";
   o << std::endl; 
   o << "Distribution functions:\n";
   o << "-l, --lhapdf=\"SET MEMBER\" : LHAPDF set and member "
     << "for the unpolarised PDF (default: GRV98)\n";
   o << "-s, --sivers=N : Sivers distribution function (default: none)\n"
     << "                 1 = Torino group 2009 fit\n"
     << "                 2 = positivity bound from unpolaried PDF\n";
   o << "-T, --transv=N : transversity distribution function (default: none)\n"
     << "                 1 = Anselmino et al., Phys. Rev. D 75, 054032 (2007)\n"
     << "                 2 = Soffer bound from GRV98 and GRSV\n"
     << "                 3 = Soffer bound from chosen unpolarised PDF and DSSV\n";
   o << "--quark_transv_soffer : equivalent to --transv=3\n";
   o << "-f, --ffunc=N : quark and gluon fragmentation type\n";
   o << "                1 = DSS\n";
   o << "                2 = Fragment parton into a single hadron\n";
   o << "      (default) 4 = DSS + Collins (pion) Torino group fit\n";
   o << "                5 = DSS + Collins (pion), positivity bound\n";
   o << "-r, --order=N : order\n";
   o << "                0 = no initial/final-state radiation,"
     << " leading order fragmentation functions\n"
     << "      (default) 1 = initial/final-state radiation,"
     << " NLO fragmentation functions\n";
   o << "-h, --hadron=N : (with ffunc=5) choose hadrons for collins fragmentation\n"
     << "               Can be repeated for multiple species e.g. -h 211 -h 111\n";
   o << std::endl; 
   o << "Physics parameters:\n";
   o << "--min_hard_pt, -m : minimum pT in pythia hard scattering "
     << "(default = " << PYTHIA_MIN_PARTON_PT << " GeV/c)\n";
   o << "-E, --eta=\"MIN MAX\" : set pseudorapidity range for hadrons\n";
   o << "--min_hadron_pt=PT, -t : minimum hadron pT (default 0)\n";
   o << "--primordial_kt=KT, -k : primordial kT of partons in hadrons "
     << "(default 0.5 GeV/c)\n";
   o << std::endl; 
   o << "Event generation (pick one):\n";
   o << "--qcd2x2 : generate QCD 2x2 scattering (default)\n";
   o << "--jpsi : J/Psi production\n";
   o << "-p, --pytune=PYTUNE : select PYTHIA tune\n";
   o << std::endl; 
   o << "Output:\n";
   o << "-e, --events=NEVENTS : number of scattering events to generate\n";
   o << "-o, --output=FILENAME : root histogram and output file\n";
   o << "-q, --qa=FILENAME : generate QA figures in a ROOT file\n";
   o << "--help : print this help screen\n";
   o << std::endl; 
   o << "Machine Parameters:\n";
   o << "--rhic200 : RHIC collider 200 GeV (default)\n";
   o << "--rhic500 : RHIC collider 500 GeV\n";
   o << "--fnalftp : FNAL E704 p+p fixed target at 200 GeV "
     << "(uses equivalent CMS energy)\n";
   o << "--fnalftpb : FNAL E704 pbar+p fixed target at 200 GeV "
     << "(uses equivalent CMS energy)\n";
   o << std::endl; 
   o << "-n            : partonhadron out 0 [default], Tparticle 1 or both 2" << std::endl;
   
}

void Configuration::Speak(std::ostream& o) const {
   o << "Output file: " << outfile << std::endl;
   o << "Fragmentation function type: " << fftype << std::endl;
   o << "Fragmentation order: " << order << std::endl;
   o << "Number of events: " << numEvts << std::endl;
   o << "Machine type: " << machine_type << std::endl;
   o << "PYTHIA event type: " << pyth_evt_type << std::endl;
   o << "Sivers set: " << sivers_type << std::endl;
   o << "Minimum hard pT: " << pyth_parton_pt << std::endl;
   o << "Mean hadron squared pT vs. parton: " << mean_jt2 << std::endl;
   o << "Hadron pseudorapidity range: " << etamin << " to " << etamax << std::endl;
   o << "Minimum hadron pT (lab): " << min_hadron_pt << std::endl;
   o << "output type: " << tparticleout << std::endl;
}

bool Configuration::Verify() {
   bool valid(true);
   // Check that the order is a valid range
   if(Order() < 0 or Order() > 1) {
      std::cerr << "Invalid order parameter:" <<
      "\n\tLO  = 0 (no pythia IS/FS interactions, LO FF)" <<
      "\n\tNLO = 1 (pythia IS/FS ON, NLO FF)" << std::endl;
      valid = false;
   } // if
   // Check that fragmentation functions that require a particular
   // order have that order.
   // Give warning but not error (valid still == true).
   if(Order() == LO and FragmentationFunction() == AKK) {
      std::cerr <<
      "WARNING: LO is requested but AKK FF's are only NLO. Calculation may not be self-consistent."
      << std::endl;
   } // if
   if(Order() == NLO and FragmentationFunction() == DSSC) {
      std::cout <<
      "WARNING: NLO is requested but Collins parametrization uses DSS LO. Order forced to 0 (LO)."
      << std::endl;
      order = LO;
   } // if
   return valid;
}

} // namespace tppmc
